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Abstract 

We have studied the dynamical properties of finite iV-unit FitzHugh-Nagumo 
(FN) ensembles subjected to additive and/or multiplicative noises, reformulating 
the augmented moment method (AMM) with the Fokker-Planck equation (FPE) 
method [H. Hasegawa, J. Phys. Soc. Jpn. 75, 033001 (2006)]. In the AMM, 
original 2A r -dimensional stochastic equations are transformed to eight-dimensional 
deterministic ones, and the dynamics is described in terms of averages and fluctu- 
ations of local and global variables. The stochastic bifurcation is discussed by a 
linear stability analysis of the deterministic AMM equations. The bifurcation tran- 
sition diagram of multiplicative noise is rather different from that of additive noise: 
the former has the wider oscillating region than the latter. The synchronization in 
globally coupled FN ensembles is also investigated. Results of the AMM are in good 
agreement with those of direct simulations (DSs). 
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1 INTRODUCTION 



The FitzHugh-Nagumo (FN) model [U [2] has been widely adopted as a simple model 
for a wide class of subjects not only for neural networks but also for reaction-diffusion 
chemical systems. Many studies have been made for the FN model with single elements 
[3]- [12] and globally-coupled ensembles [12]-[T7j. The FN model is usually solved by direct 
simulation (DS) or the Fokker-Planck equation (FPE) method. For iV-unit FN model, 
DS requires the computational time which grows as N 2 with increasing N. The FPE 
method leads to (2iV + l)-dimensional partial equations to be solved with appropriate 
boundary conditions. A useful semi-analytical method for stochastic equations has been 
proposed, taking account of the first and second moments of variables [5]. Recently we 
have proposed in [131 EB] the augmented moment method (AMM) based on a macroscopic 
point of view. In the AMM, we describe the properties of the stochastic ensembles in terms 
of a fairly small number of variables, paying our attention to their global behavior. For the 
iV-unit stochastic systems, each of which is described by K variables, ii'A-dimensional 
stochastic equations are transformed to iV^-dimensional deterministic equations in the 
AMM where N eq = K(K + 2) independent of N. This figure is, for example, N eq = 3 for 
the Langevin model (K = 1) and N eq = 8 for the FN model (K = 2). The AMM has been 
successfully applied to a study on the dynamics of coupled stochastic systems described by 
the Langevin, FN and Hodgkin-Huxley models subjected to additive noises with global, 
local or small- world couplings (with and without transmission delays) [TBI [T6l [18] . 

In recent years, the noise-induced phenomena such as stochastic resonance, noised- 
induced ordered state and noised-induced bifurcation have been extensively studied. In- 
teresting phenomena caused by additive and multiplicative noises have been intensively 
investigated (for a recent review, see Refs. [13 [2D]: related references therein). It has been 
realized that the properties of multiplicative noises are different from those of additive 
noises in some respects as follows. (1) Multiplicative noise induces the transition, creating 
an ordered state, while additive noise is against the ordering [21]-[29]. (2) Although the 
probability distribution in Langevin systems subjected to additive Gaussian noise follows 
the Gaussian, multiplicative Gaussian noise generally yields non-Gaussian distribution 
[3Q]-[33]. (3) The scaling relation of the effective strength: (3(N) = (3(l)/\/~N valid for 
additive noise is not applicable for multiplicative noise: a(N) ^ a(l)/\fN, where a(N) 
and P(N) denote effective strengths of multiplicative and additive noises, respectively, of 
iV-unit systems [34J. 
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In order to show the above item (3), the present author has adopted the AMM for the 
Langevin model in a recent paper [31] . The AMM was originally developed by expanding 
variables around their mean values in stochastic equations to obtain the second-order 
moments both for local and global variables [13] . To extend the applicability of the 
AMM to stochastic systems including multiplicative noises, we have reformulated it for 
the Langevin model with the use of the FPE [34, 35J. It has been pointed out [34J that 
a naive approximation of the scaling relation for multiplicative noise: a(N) = a(l)/\/~N, 
as adopted in [2H], leads to the result which is in disagreement with that of DSs. 

It is doubly difficult to study analytically the dynamical properties of stochastic sys- 
tems with finite populations. Most of analytical theories having been proposed so far are 
limited to infinite systems. Usually we solve the FPE for N = oo ensembles by using the 
mean-field and diffusion approximations to get the stationary probability distribution. 
For a study of dynamics, we have to obtain the instantaneous probability distribution 
from the partial differential equations (DEs) within the FPE, which is difficult even for 
N = oo. Recently, the time- dependent probability distribution is treated with a series 
expansion of the Hermite polynomials, with which dynamics of N = oo stochastic systems 
is expressed by the time-dependent expansion coefficients [12]. In the AMM, equations of 
motions for N eq moments which have clear physical meanings may describe the dynamics 
of stochastic systems with finite N. 

In this paper, we will study effects of additive and/or multiplicative noises on the 
dynamical properties of the FN model. Although effects of additive noise on the FN 
model have been extensively investigated [3] -[17], there have been no such studies on 
the effect of multiplicative noise, as far as the author is concerned. We are interested 
in the stochastic bifurcation, which is one of interesting phenomena induced by noise 
(Refs. [3Tfl I3T]. related references therein). The theory on stochastic bifurcation is still 
in its infancy. Indeed, there is no stringent definition of the stochastic bifurcation. At 
the moment, two kinds of definitions have been proposed: (i) one is based on a sudden 
change in the stationary probability distribution, and (ii) the other is based on a sudden 
change in the sign of the largest Lyapunov index. Unfortunately these two definitions 
do not necessarily yield the same result. The bifurcation of the single [UJ and ensemble 
FN model [12] [T7J subjected to additive noise has been recently discussed. Based on the 
second-order moment method, the bifurcation analysis has been made for globally-coupled 
FN model in Ref. JT7J, where dynamics of fast variables is separated from and projected 
to that of slow variables subjected to additive noise. We will discuss the bifurcation in the 
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FN ensembles subjected to additive and/or multiplicative noises, making a linear stability 
analysis to our deterministic AMM equations. It is much easier to study the deterministic 
DEs than stochastic DEs. 

The purpose of the present paper is two folds: (i) to reformulate AMM for the FN 
model subjected to additive and multiplicative noises with the use of FPE [MJ ES] , and 
(ii) to discuss the respective roles of additive and multiplicative noises on the stochastic 
bifurcation and synchronization. The paper is organized as follows. In Sec. II, we will 
apply the AMM to finite A-unit FN ensembles subjected to additive and multiplicative 
noises. With the use of the AMM equations, the stochastic bifurcation is discussed in 
Sec. III. Some discussions on the synchronization are presented in Sec. IV. The final Sec. 
V is devoted to conclusion. 



2 FN neuron ensembles 
2.1 Description of our model 

We have adopted A-unit stochastic systems described by the FN model subjected to 
additive and multiplicative noises. Dynamics of the coupled ensemble is expressed by 
nonlinear DEs given by 

^ = F(x,)-c^ + «G(x^(t)+/3^(t) + 4 c) (t) + / (e) W, (1) 

% = bxi-dyi + e, (z = ltoA) (2) 

at 

with 

4 c) (*) = f Efe-^)- (3) 

In Eq. (l)-(3), F(x) = a 3 x 3 + a 2 x 2 + a\X, a 3 = —0.5, a 2 = 0.55, a x = —0.05, b = 0.015, 
c = 1.0, d = 0.003 and e = [5J [13]: x^ and yi denote the fast and slow variables, 
respectively: G(x) is an arbitrary function of x: I^ e \t) stands for an external input: J 
expresses the strength of diffusive couplings, Z — N — 1: a and (3 denote magnitudes 
of multiplicative and additive noises, respectively, and rji(t) and express zero-mean 
Gaussian white noises with correlations given by 

(vM vAt')) = ( 4 ) 

= M(*"0, (5) 

(»*(*) 6(0) = o. (6) 
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The Fokker-Planck equation p({xi}, {yi}, t) is expressed in the Stratonovich represen- 
tation by|33] [38] 

Q Q Q 

-Q t V = ~l£^{l F ( x k)-cy k + I k ]p}-Y,Q^l( bx k- d yk + e)p\ 

+ (T)?^ G ^^ [GWri} + ?(f)S R <7) 

where I k = I [ k c) + 1^ . 

We are interested also in dynamics of global variables X(t) and Y(t) defined by 

x{t) = ^Y:x i (t), (8) 

i 

n*) = i !>(*)• ( 9 ) 

The probability of P(X, Y, t) is expressed in terms of p({xi}, {y^}, t) by 

P(X,F,t) = 1 1 n^dyi p i ({x i }dVi},t)s(x-^'Exi\ S\Y-^J2y^j 
Moments of local and global variables are expressed by 

(4Vi) = =] j dxidy y Pifay^t) x\y\, (11) 
(X k Y e ) = f f dXdYP{X,Y,t) X k Y e . (12) 



(10) 



By using Eqs. (1), (2), (7) and (11), we get equations of motions for means, variances 
and covariances of local variables by 

^ = (F{x i ))-c(y i ) + *(G r {x i )G{x i )), (13) 
d(yi) 



dt 

~dt 



dt 

d(xj yj) 
dt 



b(xi) - d(yt) + e, (14) 
(xiF(xj)) + (xjF(xi)) - c{{xiyj) + (x^)) 

+ zJ2(( X i X k) + ( X 3 X k) - ( x i) - (tf)) 
k 

a 2 

+—[(x i G'(x j )G(x j )) + (xjG'ixJGixi))] 

+ {a\G{x l f) + (15) 

Ki^Vj) + fayi)) - 2d (yiVj), (16) 

(yj F(xi)) - c{yiyj) + b(xiXj) - d(xi yj) 

9 

+^ Vj) - ( x i Vi)) + y(% G'(xi)G(xi)), (17) 



where G'(x) = dG(x)/dx. Equations (13)-(17) may be obtainable also with the use of the 
Furutsu-Novikov theorem [39, 40J. 

From Eqs. (8), (9) and (12), we get equations of motions for variances and covariances 
of global variables: 

d(V K ) 1 



dt N { 

YA = J_ TT d(v Kt vx 3 ) (Kty=12 ) (19) 

dt N 2 ~ J ~ J dt ' 1 ' 7 ' ] 1 ] 

where we adopt a convention: vu = Xi, = Hi, V\ = X and Vi = Y. Equations (13) 
and (14) were used for single (N = 1) and infinite (N = oo) FN neurons subjected only 
to additive noise (a = 0) in the mean-field approximation [12]. Equations (13)-(17) were 
employed in the moment method for a single FN neuron subjected to additive noises [5]. 
We will show that Eqs. (18) and (19) play important roles in discussing finite iV-unit FN 
ensembles. 

2.2 AMM equations 

In the AMM [13], we define eight quantities given by 

= (K)-^EW. (2°) 

T«,A = JrJ2(( V ™ -^k) {V\i -fJ-x)), (21) 

i 

p K ,x = ((V K -fi K )(V x -p x )), (M = l,2) (22) 

with 7i j2 = 72,1 and pi^ = P2,i- It is noted that j^x expresses the averaged fluctuations 
of local variables while p K ^x denotes those of global variables. 

Expanding Eqs. (13)-(19) around means of p K as v K i = pi K + 5v K i and retaining the 
terms of 0((Sv K iSvxj)), we get equations of motions for the eight quantities given by 

dn o? 

-rj- = /„ + /27i,i-cft + y[#o£i + 3 (9i92 + 9o92hi,i] + I [e \ (23) 
dfi2 



dt 



bux - d/i 2 + e, (24) 



rf 7i,i n/ __. ... N , 2JN,_ _ , , n 2 . 2 



(Jf 2(071,1 - C7i, 2 ) + — —(pi,i - 7i,i) + 2a [g l + 2^0^3)71,1 

+ a 2 gl + (3\ (25) 

^ = 2(67 1)2 -d 7a , 2 ), (26) 
dt 
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djl 2 , / A I \ ° 2 / 2 \ 

= &7l,l + i a - «)7l,2 - C72,2 + -^-(Pl,2 - 7l,2j + y (0i + 20002)71,2, (27) 

^ = 2(ap 1 , 1 -cp 1 , 2 ) + 2a 2 (0 1 2 + 20o02)pi,i + ^ + ^, (28) 

Q>2,2 

dt 



2(6pi, 2 -dp2,2), (29) 
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f/pi 2 / ,\ tt I 

= opi,i + (a - d)pi, 2 - cp 2 ,2 + y(0i + 20 O 02)pi,2, (30) 

with 

a = fi + 3/ 3 7i,i, (31) 

£ = (l/£!)FW(pO, (32) 

£ = (l/£!)G«(pi). (33) 

The original 2iV-dimensional stochastic equations given by Eqs. (l)-(3) are transformed 
to eight-dimensional deterministic equations. Equations (23)-(30) with additive noises 
only (a = 0) reduce to those obtained previously [13]. We note that in the limit of J = 0, 
AMM equations lead to 

P.,x = (M = M) (34) 

which is nothing but the central-limit theorem describing the relation between fluctuations 
in local and average variables. In the limit of N = 1, we get p K ,A = 7 k ,a, by which the 
AMM equations reduces to the five- dimensional DEs for pi, p 2 , 71,1, 72,2 and 71,2- 

Equations (l)-(3) for = a = (3 = have the stationary solution of X{ = yi = 0. 
When we assume G(x) for the multiplicative noise given by 

G(x) = x, (35) 

the AMM equations are expressed by 

Jo + /27i,i - C A*2 H — + (36) 



dt 
dp 2 

~dT 



opi - dfi 2 + e, (37) 



07ll / \ '2. J N . , 9 2 2 nl I 

— — - = 2(071,1 - C7i >2 ) H — (pi,i -71,1) + 2a 7i,i + a Pi + p , (3* 

at z 

072,2 

dt 



2(671,2-^72,2), (39) 



—7— = 071,1 + (a - «)7i,2 - C7 2 ,2 + -^-(Pi,2 - 71,2) H — , (40) 

at Z 2 



^ = 2(a ai -cp 1 , 2 ) + 2a 2 Plll + ^ + §, (41) 
^ = 2(bp lj2 -dp 2>2 ), (42) 

—j— = opi.i + (« - «)Pi.2 - cp 2 .2 H z — • (43) 

at 2 

The stochastic bifurcation of the AMM equations given by Eqs. (36)- (43) will be investi- 
gated in Sec. 3. AMM equations for a more general form of G(x) = x \ x \ s ~ l (s > 0) are 
presented in Appendix A. Some numerical results for various s values will be discussed in 
Sec. 4 (Fig. 13). 

2.3 Properties of the AMM 

Contributions from multiplicative noise have the more complicated N dependence than 
those from additive noise. Comparing the f3 2 term in d^i^/dt of Eq. (25) to that in 
dpx^/dt of Eq. (28), we note that the effective strength of additive noise of the iV-unit 
system, j3(N), is scaled by 

(3(N) = ffi (44) 



In contrast, a comparison between the a 2 terms in Eq. (25) and (28) yield the two kinds 
of scalings: 

a(N) = for p, x term, (45) 



a(N) = a(l), for 7 Xil and pi^ terms, (46) 

The relations given by Eqs. (45) and (46) hold also for d^i^/dt and dpi 2 /dt given by 
Eqs. (27) and (30). Thus the scaling behavior of the effective strength of multiplicative 
noise is quite different from that of additive noise, as previously pointed out for Langevin 
model [33]. If the relations: a(N) = a(l)/y/N and /3(N) = P(l)/\/N hold, the FPE for 
P(X, Y, t) of the global variables of X and Y in iV-unit systems may be expressed by 

^P(X,Y,t) = ~{[F(X)-cY + I}P(X,Y,t)}-^[{bX-dY+e)P(X,Y,t)} 

Unfortunately it is not the shown in Eq. (46), although Eq. (47) may be valid in 

the case of additive noise only (a = 0) [12]. 
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3 STOCHASTIC BIFURCATION 



3.1 Additive versus multiplicative noises 
3.1.1 The case of iV = 1 

In order to get an insight to the AMM, we first examine the case of single element (N — 1), 
for which AMM equations for /ix, /i 2 , 7i,i, 72,2 and 7 12 are given by 

^ = fo + h^-c, 2 + ^ + lV, (48) 

dfx 2 
~dt 

^1,1 _ ™, \ 1 0„2„, , 2 .2 1 o2 



^72,2 



6/xi — cfy^ + e, (49) 
2(a 7 i,i - c 7 i, 2 ) + 2a 2 7ljl + a 2 /i 2 + /3 2 , (50) 



2(671,2-^72,2), (51) 



= 671,1 + (a - d)7i,2 - ^2,2 + (52) 

where a = /1 + 3/371^. We have applied a step input given by 

/( £ )(t) = A9(i-t m ), (53) 

where A = 0.1, t in = 50 and 0(x) denotes the Heaviside function: O(rr) = 1 for x > 
and zero otherwise. Equations (48)-(52) for iV = 1 have been solved by using the fourth- 
order Runge-Kutta method with a time step of 0.01 and with zero initial data. Direct 
simulations (DSs) for the iV-unit FN model given by Eqs. (l)-(3) have been performed 
by using the Heun method with a time step of 0.003 and with the initial data of Xi(0) 
and 2/i(0) which are randomly chosen from [—0.01, 0.01]. Results of DS are averaged over 
1000 trials otherwise noticed. All quantities are dimensionless. 

Figure 1 shows time courses of fii(t), ^(t), 7i,i(t), 72,2(0 and 71,2(0 calculated by 
AMM (solid curves) and DS (dashed curves) with a = 0.01, (3 = 0.0 (N = 1). A single 
FN neuron fires when the external input I^(t) is applied for t > 50. By an applied 
input, Hi, H2-, 71,1, 72,2 and 7^2 show the time-dependent behavior, and they approach 
the stationary values at t > 300. Results calculated by AMM are in good agreement with 
those of DS. 

Time courses of fXi(t) [71,1 (£)] with a = 0.05, 0.1, 0.2 and 0.5 are plotted in Fig. 2(a)- 
(d) [Fig. 2(e)-2(h)] when an applied input given by Eq. (53) is applied. With increasing 
a, the magnitude of 71,1 (t) is much increased, although the profile of is almost the 
same except for a = 0.5. The time course of H\(t) for a = 0.5 in the AMM shows an 
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oscillation which is expected to be due to the stochastic bifurcation. Although the result 
of DS averaged over 1000 trials shows no oscillation, that of a single trial clearly shows 
the oscillation. 

As was shown in Fig. 2(d) and 2(h), the bifurcation may be induced by strong noise. 
In order to discuss the bifurcation, we have applied a constant input given by 

I {e \t) = I. (54) 

The stationary equations given by Eqs. (48)-(52) with d/ii/dt = et al. are solved by the 
Newton- Raphson method. Then dynamics given by Eqs. (48)- (52) is solved with initial 
values of the stationary solutions. Figures 3(a) and 3(b) show the time courses of Hi(t) 
for inputs of I = 0.1 and / = 0.5, respectively, with multiplicative noise (a = 0.1 and 
(3 = 0.0): solid and dashed curve express the results of the AMM and DS with a single 
trail, respectively. We note that for I = 0.5, /ii begins to oscillate and its magnitude is 
gradually increased while for / = 0.1, Hi shows no time development. This is more clearly 
seen in the /Ji-// 2 plots shown in Fig. 3(b) and 3(d). The oscillation for / = 0.5 is due to 
the stochastic bifurcation. Figures 3(e) and 3(f) will be explained later in Sec. 3.1.2. 

We have calculated the bifurcation transition diagrams by making a linear stability 
analysis to the deterministic AMM equations. The 5x5 Jacobian matrix C of the AMM 
equations (39)-(43) is expressed with a basis of (/xi, /i 2 , 7i,i, 72,2, 71,2) by 



C 



fo + /27i,i + « 2 /2 


— c 


h 








b 


-d 











2[(/J + 3/hi,i)7i,i + afyi] 





2(/i + 6/371,1 + « 2 ) 





-2c 











-2d 


2b 


(tt + 3/^71,1)71,2 





b + 3/371,2 


— c 


a — d + a 



,(55) 



where a = j\ + 3/371,1 and f' k = df^jd[i\. The instability is realized when any of real 
parts of five eigenvalues in the Jacobian matrix C is positive. For deterministic FN 
neuron without noises (a — (5 — 0.0), the critical condition is given by f — d = 3 + 
2 a2/ii + a\ — d = for the stationary fii, from which the oscillating state is realized for 
0.26 < / < 3.34. 



Multiplicative noise 

In order to obtain the transition diagram, we have performed calculations of the sta- 
tionary state and eigenvalues of its Jacobian matrix C, by sequentially changing a model 
parameter such as a, ft, J and /. A calculation of the stationary state for a given a 
value, for example, has been made by the Newton-Raphon method with initial values 
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which are given from a calculation for the preceding a value. Figure 4(a) shows the I -a 
transition diagram obtained for multiplicative noise (/3 = 0.0, N = 1). When changing a 
parameter as mentioned above, we have the continuous (second-order) and discontinuous 
(first-order) transitions. Solid curves denote the boundaries of the second-order transition 
between the oscillating (OSC) and non-oscillating (NONOSC). As for the first-order tran- 
sition, at / = 2.0, for example, the OSC— ^NONOSC transition takes place at a = 0.11 
when a is increased from below, while the NONOSC^OSC transition occurs at a — 0.04 
when a is decreased from above. The state for 0.04 < a < 0.11 with hysteresis is referred 
to as the OSC' state hereafter. In the OSC' state, we have two stationary solutions, as 
will be discussed below. 

Dashed curves in Fig. 5(a) show the / dependence of the maximum real part among 
five eigenvalues, \ m (referred to as the maximum index), for multiplicative noise (a = 0.1, 
P = 0.0). Two dashed curves for 0.19 < I < 2.29 express the results of the two stationary 
solutions in the OSC state shown in Fig. 4(a). The lower, dashed curve crosses the zero 
line at four points at I = 0.29, 1.41, 2.39 and 3.41. For a comparison, the result for the 
deterministic model (a = (3 = 0.0) is plotted by the chain curve. Although the OSC state 
disappears for fairly strong multiplicative noises of a > 0.16 at 2 ~ / ~ 3.3, it persists to 
strong noises at ~ / ~ 1.3. 

Solid and dashed curves in Fig. 6(a) show fj,i and 7^1, respectively, of the stationary 
state as a function of I. Two dashed curves for 0.19 < I < 2.29 express 71,1 of the two 
stationary solutions in the OSC' state. In contrast, there is little difference in fii for 
the two stationary solutions. We note that /ix is nearly proportional to / and that the 
upper curve of 7^1 has a broad peak centered at / ~ 2. Open and filled circles in Fig. 
4(a) denote two sets of parameters of {I, a) = (0.1,0.1) and (0.5,0.1) adopted for the 
calculations shown in Figs. 3(a) and (b). 

Additive noise 

Figure 4(b) shows the I-f3 transition diagram for additive noise (a = 0.0, N = 1). 
Solid curve expresses the boundary between the OSC and NONOSC states, and dashed 
curve denotes the boundary of the first-order transition relevant to the OSC' state with 
hysteresis. When an additive noise is included, the boundary for the OSC-NONOSC 
states is much modified. This is explained for the case of /3 — 0.1 in Fig. 5(b), where the 
chain curve denotes the maximum index A m for the case of a = f3 = 0.0 while the dashed 
curve shows the result for additive noise of (5 — 0.1 (a = 0.0). The chain curve crosses 
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the zero line at I = 0.26 and 3.34, while the dashed curve has the zeros at four points 
at / = 0.12, 0.86, 2.75 and 3.48. The transition diagram has a strange shape, which is 
symmetric with respect to the axis of / = 1.80. The / dependences of jii and 7^1 are 
plotted by solid and dashed curves, respectively, in Fig. 6(b). 7^1 for additive noise has 
a broad peak similar to that for multiplicative noise shown in Fig. 6(a). It is noted that 
•y/71,1 expresses the effective width of the probability distribution of p(x) = J p(x, y) dy. 
We note that 7^1 in Fig. 6(a) shows a rapid increase at I ~ 0.2 where the OSC-NONOSC 
transition takes place in Fig. 4(a). Except this case, however, there are no abrupt changes 
in 7i t i at the transition. 

3.1.2 The case of finite N 

Multiplicative noise 

Figure 3(e) shows the time courses of /xi(t) for an input of / = 0.5 given by Eq. 
(54) applied to N = 100 FN ensembles with J = 1.0 subjected to multiplicative noise 
(a = 0.1, /3 = 0.0): solid and dashed curves express the results of the AMM and DS 
with a single trail, respectively. It shows that fii begins to oscillate and its magnitude 
is rapidly increased. This is more clearly seen in the 1^1-/^2 plots shown in Fig. 3(f). In 
contrast, for smaller J = 0.1, we have not obtained the oscillating solution, just as in the 
N — 1 case shown in Figs. 3(a) and 3(b). 

We have performed a linear analysis for the case of finite N, by using the 8x8 Jacobian 
matrix with the stationary solutions obtained from Eqs. (36)-(43), expressions for 8 x 8 
Jacobian matrix elements being presented in Appendix B. 

The I-ot transition diagram for multiplicative noise with N = 100 and J = 1.0 is 
shown in Fig. 4(c). The region of the OSC state is a little decreased for a weak a but it 
is widen for stronger a. The solid curve in Fig. 5(a) expresses the / dependence of the 
maximum index A m for multiplicative noise with a = 0.1 and J = 1.0, which crosses the 
zero line at two points at / = 0.21 and 3.37. 

Figure 7(a) shows the J-a transition diagram for multiplicative noise {(3 = 0.0) for 
/ = 3.0 with iV = 100. For a = 0.1, the OSC state is realized even for J = 0.0. For strong 
multiplicative noise of a = 0.3, the OSC state disappears at J < 0.365. For a = 0.2, the 
OSC state is realized not only at J > 0.194 but also at 0.085 < J < 0.136, as shown in 
the inset of Fig. 7(a). This implies the re-entrance to the OSC state from the NONOSC 
state when the coupling is decreased from above. 
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Additive noise 

Figure 4(d) shows the I- (3 transition diagram for additive noise for J = 1.0 and 
iV = 100. The width of the OSC state is gradually decreased with increasing (3. Chain 
and solid curves in Fig. 5(b) show the / dependences of A m for (3 = 0.0 and (3 = 0.1, 
respectively, with J = 1.0 and N = 100. The former has the zeros at / = 0.26 and 3.34 
while the latter at / = 0.29 and 3.32. 

Figure 7(b) shows the J-f3 transition diagram for additive noise for / = 3.0 with 
N = 100. The critical values of (3 are 0.114, 0.221 and 0.265 for J = 0.0, 0.5 and 1.0, 
respectively. The region of the OSC state is increased with increasing J. 

When we compare the bifurcation transition diagrams for multiplicative noise in Fig. 
4(a) and 4(c) with those for additive noise in Fig. 4(b) and 4(d), we note that the former 
is rather different from the latter, in particular for the N — 1 case. When a weak additive 
noise is added to the N — 1 model, the OSC state is slightly increased although for a large 
noise, the OSC state disappears. In contrast, the OSC state persists for multiplicative 
noise. Figures 7(a) and 7(b) show that the coupling is beneficial to the OSC state for 
both additive and multiplicative noises, as expected. 

3.2 Coexistence of additive and multiplicative noises 

Although we have so far discussed the additive and multiplicative noises separately, we 
now consider the case where both the noises coexist. Figure 8(a) shows the I — a transition 
diagram for (3 = 0.05 and N — 1. It is similar to the transition diagram shown in Fig. 
4(a) for the case of multiplicative noise only ((3 = 0.0). In Fig. 8(a), the OSC state is 
completely split with a gap even for a = 0.0. This is because the OSC state disappears 
for additive noise of (3 = 0.05 even without multiplicative noise as shown in Fig. 4(b). 

In contrast, Fig. 8(b) shows the I — 13 transition diagram for a — 0.1 and N = 1. We 
have the OSC state for weak (3 at 1.30 < I < 2.38. This is related to the fact that the 
OSC state exists for a = 0.1 and (3 = 0.0 in Fig. 4(a). 

Figure 9 shows the a-f3 transition diagram for the OSC and NONOSC states with 
/ = 1.0 and / = 3.0 (N = 1). In the case of / = 1.0, the boundary between the OSC and 
NONOSC states extends to large a, reflecting the behavior shown in Figs. 4(a) and 8(a). 
We have tried to fit the calculated boundaries by a simple expression of (3 = c^Jl — (a/d) 2 . 
Dashed curves in Fig. 9 express the results with c = 0.084 and d = 0.32 for / = 1.0 and 
with c = 0.121 and d = 0.158 for / = 3.0. An agreement between the solid and dashed 
curves in the case of / = 3.0 is better than that in the case of / = 1.0. 
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4 DISCUSSION 



It is interesting to study the synchronization in FN ensembles with noises. In order to 
quantitatively discuss the synchronization in the ensemble, we first consider the quantity 
given by [13] 

R (t) = 4 E<M*) - *i(*)] 3 > = 2 bi,iW - (56) 

When all neurons are in the completely synchronous state, we get Xi(t) = X(t) for all 
i, and then R(t) = in Eq. (56). On the contrary, in the asynchronous state, we get 
R(t) = 2(1 - 1/A07M = ^o(*) because p lyl = j hl /N [Eq. (34)]. We modify R(t) such 
that the synchronization is scaled between the zero and unity, defining the synchronization 
ratio given by 



Ro(t) 



N-l 



(57) 



which is and 1 for completely asynchronous (R = Rq) and synchronous states (R = 0), 
respectively. As will be shown below, S(t) depends not only on model parameters such 
as J and iV but also the type of noises (a and (5). 

The calculations have been performed for an external input I^ e \t) given by Eq. (53). 
DS calculations have been made with 20 trials. 

Figures 10 (a)- (d) show time courses of fi\(t), 7i,i(t), Pi,i(£) and S(t) for multiplicative 
noise (a = 0.01, f3 = 0.0) with J = 1.0 and iV = 100. In contrast, Figs. 10(e)-(h) express 
the result for additive noise (a = 0.01, (3 = 0.0). When the external input is applied 
at t — 50, FN neurons fire, and 7i,i(£), Pi,i(t) and S(t) develop. The time dependence 
of Hi(t) for multiplicative noise in Fig. 10(a) is almost the same as that for additive 
noises shown in Fig. 10(e). From a comparison between Figs. 10(b) and 10(f), we note, 
however, that 7^1 for multiplicative noise is considerably different from that for additive 
noise. This is true also for px,i shown in Figs. 10(c) and 9(g). The time course of S(t) 
is calculated with the use of Eq. (57) with 711 and pi t i shown in Figs. 10(b) and 10(c) 
[Figs. 10(f) and 10(g)]. Reflecting the differences in 7^ and p^i, the result of S(t) for 
multiplicative noise in Fig. 10(d) is quite different from that for additive noise in Fig. 
10(h), although at t > 400, both S(t) approach the same value of S = 0.24. 

We may apply also spike train and sinusoidal inputs to the systems, whose results will 
be discussed in the followings. Figures 11(a) and 11(b) show time courses of pi(t) and 
S(t) in the case of multiplicative noise (a = 0.01, (3 = 0.0) with J = 1.0 and N = 100, 
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when a spike train input given by 

I^ e \t) = A^(t-t n )e(t n + T w -t), (58) 

n 

with A — 0.1, T n — 50 + 100 (n — 1) and T w = 10 is applied, an input being plotted at 
the bottom of Fig. 11(a). In order to understand the relation between ^ and S, we 
depict the S — (dfxi/dt) plot shown by solid curves in Fig. 11(c), where initial data of S 
and d/ii/dt at t < 350 are neglected. We note that S is large for d/ii/dt ~ —0.15 and 
< d/Ji/dt < 0.1. In contrast, we show by the dashed curve in Fig. 11(c), a similar plot 
in the case of additive noise (a = 0.0, (3 = 0.01). The trend of multiplicative noise is 
similar to that of additive noise for d/ii/dt < but quite different for d/ii/dt > 0. 

Figures 12(a) and 12(b) show time courses of /xi(t) and S(t) in the case of multiplicative 
noise (a = 0.01, (3 = 0.0 with J = 1.0 and N = 100, when a sinusoidal input given by 

'2n(t-t b y 



&\t) = AQ(t-t b ) 



cos 



T 



(59) 



with A — 0.1, t b — 50 and T p = 100 is applied, input being plotted at the bottom of Fig. 
12(a). The solid curve of Fig. 12(c) expresses the S — (dfxi/dt) plot, where initial data of 
S and dfii/dt at t < 350 are neglected. In contrast, the dashed curve in Fig. 12(c) shows 
the S — (d/ii/dt) plot in the case of additive noise (a = 0.0, (3 = 0.01). It is noted that 
the behavior of the S — (d/ii/dt) plot of multiplicative noise is similar to that of additive 
noises for d/ii/dt < but rather different for d/ii/dt >. This is consistent with the result 
shown in Fig. 11(c) for spike train. It is interesting that the synchronization S seems to 
correlate with | d\i\jdt |: S becomes larger for a larger | d/ii/dt \. 

We have discussed the stochastic bifurcation and the synchronization in FN ensembles 
subjected to additive and/or multiplicative noises. Our calculations have shown that 
effects of multiplicative noise are rather different from those of additive noise. This may 
be understood by a simple analysis as follows. Equation (48) is rewritten as 

d\x o? 

— - = a 3 nl + a 2 f4 + ai^i - cfi 2 + (3a 3 //i + a 2 )7i,i + yAn + I- (60) 

For weak noises, we get 

7i,i ~D(a 2 fii + (3 2 ), (61) 

where D is the coefficient to be determined from Eqs. (48)- (52) but its explicit form is 
not necessary for our discussion. Substituting Eq. (61) to Eq. (60), we get 

— a' 3 nl + a' 2 fx 2 + a'-tHi - cfi 2 + I' ', (62) 



dt 
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with 



a. 



I 



a 3 (l + 3Da 2 ), 
a 2 (l + Da 2 ), 

2 

a x + y + 3Da 3 p 2 , 
I + Da 2 (3 2 . 



(63) 
(64) 

(65) 
(66) 



The dynamics of fii and \i 2 is effectively determined by Eqs. (49) and (62). Results 
for deterministic FN model are given by setting a — (3 — 0.0 in Eqs. (62)-(66). We 
note in Eqs. (63)-(66) that additive noise modifies constant and linear terms, whereas 
multiplicative noise changes the linear, quadratic and cubic terms. These differences yield 
the difference in the stochastic bifurcations for additive and multiplicative noises. 

The bifurcation diagram of a single FN model with additive noise is discussed in 
Refs. [SI [HI H2J- It has been shown that the probability distribution of single FN model 
obeys Gaussian for weak additive noise while for strong noise, it shows a deviation from 
Gaussian with the bimodal structure (TTj [T5] . The transition from the Gaussian to non- 
Gaussian distribution is considered to show the stochastic bifurcation. However, a change 
in the form of probability density is generally gradual when a parameter of the model is 
changed. By employing the second-order moment method, Tanabe and Pakdaman [11] 
have obtained the bifurcation diagram showing a critical current I c of a single FN model 
with additive noise, which approaches I c for the deterministic FN model with decreasing 
the strength of additive noise. Our transition diagram shown in Fig. 4(b) is different 
from their transition diagram (Fig. 4 of Ref. |llj ) in which the OSC state exists even for 
a strong additive noise. By solving the FPE by an expansion of the Hermite polynomials, 
Acebron, Bulsara and Rappel [12] have shown that the stochastic bifurcation cannot be 
obtained for a single FN model against our result showing the bifurcation for N — 1 [Fig. 
4(b)]. Recently, stochastic bifurcations in globally coupled (N = oo) ensembles subjected 
to additive noise have been discussed in Refs. [121 EE Ej- It nas been shown [12] 
that contrary to a single element, the bifurcation occurs in globally coupled ensembles 
as the noise strength is increased. With the use of the moment method, Zaks, Sailer, 
Schimansky-Geier and Neiman [17] have discussed the bifurcation of coupled (N = oo) 
FN model in which slow variables [yi) are subjected to additive noise whereas in our 
study, fast variables (xj) are subjected to additive and/or multiplicative noises [Eqs. (1) 
and (2)]. With increasing the noise intensity, mean field shows a transition from a steady 
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equilibrium to global oscillations, and then back to another equilibrium for sufficiently 
strong noise [TT] . 

In the conventional moment (or cumulant) approach, the Gaussian distribution is as- 
sumed for calculations of first and second moments, and the moment method is considered 
to lose its validity for the non-Gaussian distribution [5]-[II],[I7]. Our reformulation of 
the AMM with FPE has revealed that the moment method is free from the Gaussian 
approximation and that it is valid also for non-Gaussian distribution. Indeed, we have 
shown in Refs. [3H, [35] that the AMM can be well applied to the Langevin model with 
multiplicative Gaussian noise although its probability distribution generally follows the 
non-Gaussian [301 [31]. The moment approach is expected to have the wider applicability 
than having been considered so far. 

The transition diagrams shown in Figs. 4(a), 4(c) and 4(e) for multiplicative noise 
have the asymmetric, peculiar structure, which arises from the assumed form of G(x) = x 
in Eq. (35). If we alternatively assume G(x) = 1, 'multiplicative noise' reduces to 
additive noise, and its transition diagrams are given by Figs. 4(b), 4(d) and 4(f). Thus 
the structure of the bifurcation transition diagram depends on the adopted form of G(x). 
Some preliminary calculations have been made with the use of a form of G(x) = x \ x 
by changing the index s: the relevant AMM equations are presented in Appendix A. 

Figure 13(a) shows time courses of S(t) for various s with a = 0.01, (3 = 0.0 and 
N = 100. At a glance, the overall behavior of S(t) seems almost independent of s. We 
note, however, that 7^1 and pi i, which are plotted in Figs. 13(b) and 13(c), respectively, 
much depend on the s value. It has been shown that in the Langevin model, the stationary 
distribution shows much variety depending on the index s in G(x) = x \ x [34], 
[35] . This is expected to be true also for the distribution p(x,y) in the FN model with 
multiplicative noises. It would be interesting to investigate the stochastic bifurcation by 
changing the index s, whose study is under consideration. 

5 CONCLUSION 

We have studied effects of additive and multiplicative noises in single elements and 
globally-coupled ensembles described by the FN model, by employing AMM reformu- 
lated with the use of FPE [MJ [13]. The stochastic bifurcation has been examined by a 
linear analysis of the deterministic AMM equations. The properties of the multiplicative 
noise in FN neuron ensembles is summarized as follows. 
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(a) the effect of multiplicative noise on the stochastic bifurcation is different from that of 
additive noise, and 

(b) the effect of multiplicative noise on the synchronization is more significant than that 
of additive noise. 

The item (a) is consistent with the results obtained for other nonlinear systems such as 
Duffing- Van der Pol model J37J HI]- The effect of noise depends on the type of noises, 
and it is also model dependent |42j. The item (b) is similar to the item (1) of an ordered 
state created by multiplicative noises [21]- [29] mentioned in the introduction. 

A disadvantage of our AMM is that its applicability is limited to weak-noise cases. In 
physics, there are many approximate methods which are valid in some limits but which 
provide us with clear physical picture beyond these limits. One of such examples is 
the random-phase approximation (RPA) which has been widely employed in solid-state 
physics. The RPA is valid in the limit of weak interactions. The RPA is, however, used 
for larger interactions leading to the divergence in the response function, which signifies 
the occurrence of excitations such as spin waves. We expect that the AMM may be such 
an approximate method, yielding meaningful qualitative result even for strong noises. 

On the contrary, an advantage of the AMM is that we can easily discuss dynamical 
properties of the finite iV-unit stochastic systems. For iV-unit FN neuronal ensembles, 
the AMM yields the eight-dimensional ordinary DEs, while DS and FPE require the 
2iV-dimensional stochastic DEs and the (2N + l)-dimensional partial DEs, respectively. 
Furthermore the calculation of the AMM is much faster than DSs: for example, it is 
about 2000 times faster than DS calculations with 100 trials for 100-unit FN ensembles. 
We hope that the AMM may be applied to a wide class of coupled stochastic systems 
subjected to additive and/or multiplicative noises. 
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Appendix A: AMM equations for G(x) 

In the case of G(x) given by 



G(x) = x | x \ s - 1 {s > 0), (Al) 
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AMM equations (51)-(58) become 



dt 

dfi 2 
~dt 



fo + /27i,i - C P2 + ~a 2 s/j,i | ^ \ 2s 2 [1 + (s - l)(2s - 1) | fi x ] 2 7 M ] 
J (e) , (A2) 
- rf/i 2 + e, (A3) 

2(071,1 - C7i, 2 ) + - 7i,i) + 2a 2 s(2s - 1) | fix \ 2s ' 2 71,1 

+ a 2 \fix \ 2s +P 2 , (A4) 
= 2(671,2-^72,2), (A5) 



<*7i,i _ N , 2JA r .. > .2 , . r, 12,-2 

dt 

^72,2 

= 071,1 + (a - 0)71,2 - C72,2 + -^-(Pi,2 - 7l,2) 

+ ia 2 s(2s - 1) I fix \ 2s ~ 2 71,2, (A6) 

^ = 2(ap 1 , 1 -cp 1 , 2 ) + 2a 2 S (2 S -l)^ 2 p 1 , 1 + ^^ + ^, (A7) 

dt 

l, l>\:2 , , ,1 -1 . 1 \ I , , 2,-2 



2(6pi,2-rfp 2 , 2 ), (A8) 



= 6/01,1 + (a - d)px,2 - cp 2 ,2 + -a 2 s(2s - 1) | fix \ 2s p 1:2 - (A9) 

Some numerical examples of 711 (i) et al. for various values of the index s are shown in 
Fig. 13. 



Appendix B: Jacobian matrix of AMM equations 

For an analysis of the stochastic bifurcation of finite iV-unit FN ensembles, we need 
the 8x8 Jacobian matrix C expressed with a basis of (fix, (J>2, 7i,i, 72,2, 71,2, Pi,i, P2,2, 
Pi,2), by 



ft + #yi,i + « 2 /2 


— c 


h 





6 


-d 








2[(/i + 3/hi,i)7i,i + « 2 pi] 





2(/i + 6/371,1 + « 2 ) - 2 JN/Z 














-2d 


(fi + 3/^71,1)71,2 





b + 3/371,2 


—c 


2(/{ + 3^ 7 i,iKi + 2aWiV 





6/sPi,i 

















(/{ + 3/^l,l)pi,2 





3/ 3 Pl,2 
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-2c 2JN/Z 

2b 

a-d + a 2 /2- JN/Z JN/Z 

2(a + a 2 ) -2c 

-2d 2b 

b -c a-d + a 2 /2 



(Bl) 



where a = f\ + 3/371,1 and = dfk/dfii, the (l-4)-column and (5-8)-column elements 
being separately expressed for a convenience of printing. 
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